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Abstract 

It  has  been  shown  that  a  three-point  second  difference  estimator  is  nearly  optimal  for  estimating 
frequency  drift  in  many  common  atomic  oscillators.  We  derive  a  formula  for  the  uncertainty  of  this 
estimate  as  a  function  of  the  integration  time  and  of  the  Allan  variance  associated  with  this  integration 
time. 


Theory 


The  three-point  drift  estimator  is  a  useful  tool  for  estimating  the  frequency  drift  in  many  atomic 
oscillators  [1] .  In  this  paper  we  derive  a  formula  for  the  uncertainty  of  the  three  point  drift  estimate; 
as  we  shall  demonstrate,  there  is  a  simple  relationship  between  the  uncertainty  of  the  drift  estimate 
and  the  Allan  variance  of  the  residuals  which  remain  after  the  estimated  drift  is  removed.  We  explain 
how  to  apply  the  uncertainty  formula  and  then  we  use  it  to  assess  the  uncertainty  of  the  drift  estimate 
in  several  examples. 


Let  us  begin  by  discussing  the  three-point  drift  estimator.  To  define  it,  let  x(t)  be  a  time  series  of 
time  difference  measurements  between  two  oscillators  drifting  in  frequency  relative  to  each  other.  An 
optimal  estimator,  D,  of  drift  uses  the  first,  middle,  and  last  time-difference  points.  We  estimate  the 
average  frequency  over  the  first  and  second  halves  of  the  data,  subtract  the  first  frequency  from  the 
second,  and  then  divide  by  r,  the  time  elapsed  between  the  first  and  middle  or  middle  and  last  data 
points.  This  yields: 

-  _  1  /  x(2  t)  —  x(r)  x(t)  —  z:(0) 

T  \  T  T 

(1) 


=  ((«(2r)  -  2*(r)  +  *(<>)) 


That  is,  we  estimate  drift  as  1  /r2  times  the  second  difference  of  the  time  series  x,  where  we  take  the 
second  difference  over  as  large  an  interval  as  possible. 
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Let  us  separate  the  time  offset  x{t)  into  the  part  due  to  the  frequency  drift  D  and  the  part  due  to 
everything  else  (initial  offsets,  stochastic  noise,  systematics) : 

x(t)  =  x'(t)jt2  .  (2) 

We  now  will  show  that  the  uncertainty  of  the  drift  estimate,  D,  is  functionally  related  to  the  Allan 
variance  of  the  x'(t)  time  series. 

If  we  substitute  (2)  into  (1)  we  obtain 

b  =  ^{x'(2T)-2x'(t)+x’(0)  +  DT2')  .  (3) 

Rearrangement  yields: 

D  -  D  =  (x'(2t)  —  2 x’(t)  +  2/(0))  .  (4) 


The  expected  variance  of  our  drift  estimate,  D ,  around  the  true  drift  D  will  thus  be 

({b  -  D)2)  =  ~((x'(2r)  -  2x’(t)  +  x’(0))2)  ,  (5) 

where  (  }  is  the  expectation  operator.  The  square  root  of  this  quantity  is  the  expected  deviation  of  D 
around  the  true  value. 

If  we  compute  an  Allan  variance  of  x1  for  the  integration  time  r  we  obtain  [2,3] 

<vM!  =  <(*W  -  2x'(t)  +  A0)f)  ■  (6) 

Substitution  of  (6)  into  (5)  yields  our  result,  the  relationship  between  the  expected  deviation  in  the 
drift  estimator,  D,  and  the  Allan  variance  of  x\  the  drift-removed  data: 


(7) 


where  ay  (r)  is  the  Allan  variance  of  the  x’  data,  and  y  refers  to  the  frequency  data  derived  from 
2'.  Thus  we  see  that  the  uncertainty  in  our  drift  estimate  is  a  function  of  the  Allan  variance  of  the 
drift-removed  data. 

Application  of  Equation  7 

The  application  of  (7)  requires  a  bit  of  finesse.  First  of  all,  the  alert  reader  has  probably  noticed 
that,  since  we  don't  know  the  value  of  D,  the  true  drift,  we  cannot  obtain  the  time  series  x‘(t).  To 
circumvent  this  problem,  we  obtain  an  approximation  of  x'(t)  by  removing  the  estimated  drift  from 
the  x(t)  series.  We  then  compute  the  Allan  variances  for  the  approximate  x'(t)  series.  However,  it  is  at 
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this  point  that  we  encounter  another  problem:  It  is  generally  true  that  if  you  1)  use  a  second-difference 
estimator  (such  as  the  three-point  estimator)  to  estimate  drift,  2)  remove  this  estimated  drift  from  the 
time  series,  and  then  3)  compute  the  Allan  variances  for  the  residual  time  series,  the  Allan  variances 
obtained  for  large  integration  times  (such  as  r  =  1/2  the  data  length)  will  be  biased  low,  i.e.  the 
Allan  variance  will  not  be  an  accurate  measure  of  the  frequency  variability  at  large  integration  times. 
In  fact,  if  we'  were  to  take  a  data  set  with  constant  drift,  compute  D  using  (1),  remove  I)t2 3 4  j‘2  from 
each  data  point,  and  then  compute  f7,y(r)  for  this  same  r,  we  would  obtain  exactly  zero. 

We  need  to  have  <Jy<{r)  for  =1/2  the  data  length  in  order  to  use  (7).  Yet  we  know  that  after 
removing  D  from  x(t),  we  are  going  to  get  the  incorrect  value  of  0  for  oy‘2{r)  for  r  =  1/2  the  data 
length.  However,  while  cr^(r)  is  incorrectly  low  for  large  r,  it  does  accurately  represent  the  frequency 
variability  for  smaller  r.  Furthermore,  the  noise  processes  of  atomic  oscillators  are  such  that,  for  a 
given  range  of  integration  times  r,  it  is  usually  the  case  that  ^y(T)  —  krn,  where  k  is  a  constant  and 
n  is  an  integer  ranging  from  1  to  —2.  The  result  of  this  power-law  behavior  of  oy(r)  is  that  log-log 
plots  of  oy(r)  versus  t  exhibit  linear  behavior.  This  can  be  seen  in  Figure  1.  Therefore,  in  order  to 
obtain  0y(r)  for  r  =  1/2  the  data  length,  we  look  at  the  log-log  plot  of  cry  (r))  versus  r  and  discard 
the  incorrectly-low  values  of  cry(r)  which  occur  at  large  r  (For  example,  in  Figure  1,  we  would  discard 
the  point  for  which  log  r  (seconds)  %  7.  In  Figure  3  we  would  discard  the  point  for  which  log  r 
(seconds)  6.75).  Then,  we  use  the  r7y)r)  points  which  correspond  to  the  largest  remaining  r  values 
to  determine  k  and  n  (i.e.,  we  determine  the  equation  of  the  line  on  the  log-log  plot  formed  by  the 
remaining  valid  data  points).  Then,  knowing  k  and  n,  we  use  the  equation  a2, (r)  =  krn  to  determine 
the  value  of  Oy(r)  at  r  =  1/2  the  data  length.  This  value  is  what  we  need  to  apply  (7). 

For  cesium  beam  and  rubidium  gas-cell  oscillators,  the  dominant  noise  types  at  large  integration 
times  are  flicker  frequency  modulation  and  random  walk  frequency  modulation  (FLFM  and  RWFM, 
respectively).  FLFM  corresponds  to  an  n  value  of  0  and  RWFM  corresponds  to  an  n  value  of  +1. 
For  very  large  r,  RWFM  generally  dominates.  Therefore,  if  the  last  (i.e.  largest  r)  valid  linear  trend 
that  we  see  on  the  log-log  plot  is  consistent  with  a  model  of  RWFM,  we  may  use  this  slope  with  a 
measure  of  confidence  to  estimate  the  value  of  at  r  =  1/2  the  data  length.  If,  however,  the  last 

linear  trend  corresponds  to  FLFM,  we  need  to  ask  ourselves  whether  the  FLFM  noise  type  continues 
out  to  r  =  1/2  the  data  length,  or  whether  RWFM  is  the  correct  noise  type  for  r  =  1/2  the  data 
length.  The  assumption  of  RWFM  as  the  noise  type  always  leads  to  a  larger  computed  value  of  a'2,  (r) 
than  the  assumption  of  FLFM.  Thus,  simply  assuming  that  RWFM  dominates  at  r  =  1/2  the  data 
length  yields  a  conservative  estimate.  The  uncertainty  in  the  Allan  variance  estimate  will  limit  the 
accuracy  of  our  uncertainty  estimate.  Nevertheless,  we  can  make  conservative  estimates  of  uncertainty 
and  obtain  meaningful  results. 

In  summary,  to  use  (7)  to  estimate  the  uncertainty  of  we  take  the  following  steps: 

1.  Compute  using  the  second  difference  estimator  (3),  where  in  that  equation,  r  =  the  time 
interval  for  one-half  the  data  length.  Remove  Dt2 /2  from  each  of  the  time-difference  data  points 
x(t). 

2.  Compute  the  Allan  deviations  oy  (r),  for  r  =  nro,  where  n  is  an  integer  multiple  of  the  sampling 
interval  To-  Make  a  log-log  plot  of  a y(r)  versus  r. 

3.  Look  for  abnormally  low  values  of  cry(r)  at  large  values  of  r.  Discard  them. 

4.  Determine  the  parameters  k  and  n  in  the  equation  a2,  (r)  =  fern  for  the  last  valid  linear  trend 
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on  the  log-log  plot.  Then  use  this  equation  to  compute  o^,(r)  for  tv^.  Remember  to  consider 
the  possibility  that  the  noise  type  might  change  past  the  last  valid  oy(r)  value  on  the  log-log 
plot  (i.e.,  the  noise  type  might  change  from  FLPM  to  RWFM). 

5.  Substitute  this  value  of  (t )  into  (T).  Solve  (7)  for  the  variance  of  D.  The  square  root  is  the 

expected  deviation. 


Examples 

As  examples  we  use  atomic  standards  aboard  GPS  satellites  studied  from  July  1,  1991,  to  September 
15,  1992,  a  period  of  443  d.  Satellites  are  referred  to  by  their  pseudo-random  code  number  (PRN), 
the  number  by  which  users  identify  satellites  ,  or  by  their  satellite  vehicle  number  (SVN),  the  number 
used  by  the  GPS  control  segment.  Clocks  on  the  GPS  satellites  are  measured  at  NIST  against  the 
ATI  time  scale.  For  clocks  which  ran  for  this  entire  period,  drift  could  be  estimated  using  a  second 
difference  with  r=221.5  d.  Not  all  clocks  analyzed  were  on  line  foT  this  entire  period,  in  which  case 
shorter  r  values  were  found.  We  found  an  assortment  of  dominant  noise  types  at  various  integration 
times,  with  FLFM,  and  RWFM  dominating  at  times  equal  to  one-half  the  data  length.  Table  I  gives 
our  example  results  and  indicates  associated  figure  numbers. 

PRN#2  and  figure  1  illustrate  the  difficulty  in  determining  noise  type.  Looking  at  figure  1,  we  see 
that,  while  FLFM,  r°,  is  the  probable  slope  for  the  last  valid  a v(r)  values,  the  uncertainty  allows  for 
the  possibility  of  a  r1/2  slope,  indicating  RWFM.  Furthermore,  RWFM  is  usually  the  dominant  noise 
process  for  cesium  frequency  standards  at  integration  times  such  as  221.5  d  [6].  We  compute  a  more 
conservative  value  in  the  second  line  of  the  table.  Similarly  for  PRN#25  we  have  assumed  FLFM  in 
its  first  line.  If  we  assume  RWFM  we  see  we  find  only  a  small  change. 

Another  consideration  is  that  equation  (7)  applies  to  the  Allan  variance,  not  the  modified  Allan 
variance.  In  figures  1  and  5  we  used  the  modified  Allan  variance.  We  can  account  for  this  as  follows. 
Asymptotically,  if  we  define 


R 


co  — 


lim  mod  <t2(t) 
t-+oo  tr2(r)  ’ 


(8) 


then  Rq o  =  0.91  for  RWFM  and  =  0.82  for  FLFM  [5].  These  corrections  have  been  included  in 
the  table. 


Conclusions 

We  have  derived  a  relationship  that  allows  us  to  estimate  the  uncertainty  of  the  three-point  estimator 
of  frequency  drift.  It  does  not  give  a  lot  of  precision  but  it  is  adequate  for  determining  a  confidence 
level.  With  the  procedure  outlined,  we  can  determine  an  upper  bound  on  the  uncertainty  of  the 
estimate  of  frequency  drift. 
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Table  1 


PRN# 

(Type) 

Dates 

^"mix 

days 

Dominant 
Noise  Type 

0y(r)  @  r(s) 

Estimated  Drift  ±  o 
parts  in  10 _15/d 

Figure 

Nos. 

2  (Cs) 

1Ju191  -  15Sep92 

221.5 

FLFM 

0.4 -1013 

-2.7  ±0.3 

1 

2  (Cs) 

IJui91  -  15Sep92 

221.5 

RWFM 

0.2-1013®  106 

-2.7  ±0.6 

1 

3  (Rb) 

ljul91  -  15Sep92 

221.5 

RWFM 

2.0  *  10 13  @  106 

-98  ±6 

2 

12  (Rb) 

8Apr  -  15Sep92 

80.5 

RWFM 

2,5 -lO13®  106 

-130  ±10 

3 

19  (Cs) 

lJul  -  18Dec91 

85.5 

RWFM 

1.2-1013®  106 

35  ±5 

4 

25  (Rb) 

30Jun  -  15Sep92 

39 

FLFM 

0.7- 1013 

-183  ±3 

5 

25  (Rb) 

30Jun  -  15Sep92 

39 

RWFM 

0.6-1013®  106 

-183  ±4 

5 

PRN  2  C SVN  13 )  -  NIST (ATI) 

I  July  'SI  -  15  Sep  '92 
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Figure  1:  The  modified  Allan  variance  of  the  Cs  clock  on  PRN//2  as  measured  at 
HIST  against  the  ATI  time  scale  from  July  1,  1991  to  September  15,  1992. 
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PRN  3  < SVN  ll>  -  N  1ST (ATI) 
1  July  '91  -  15  Sep  '92 
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Figure  2:  The  Allan  variance  of  the  Cs  clock  on  PRN# 3  as  measured  at  NIST 
against  the  ATI  time  scale  from  July  1,  1991  to  September  15,  1992. 

PRN  12  (SVN  10)  Rb  Clack  -  NIST(ATl) 
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Figure  3:  The  Allan  variance  of  the  Rb  clock  on  PRN//12  as  measured  at  NIST 
against  the  ATI  time  scale  from  April  8,  1991  to  September  15,  1992. 


PRN  19  CSVN  19)  Cj  Clock  -  NI5T ( AT  1 ) 
1  July  '91  -  10  Dec  '91 


LOG  TAU  (Seconds) 


Figure  4:  The  Allan  variance  of  the  Cs  clock  on  PRN# 19  as  measured  at  NIST 
against  the  ATI  time  scale  from  April  8,  1991  to  September  15,  1992. 
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Figure  5:  The  Allan  variance  of  the  Rb  clock  on  PRN#25  as  measured  at  NIST 
against  the  ATI  time  scale  from  April  8,  1991  to  September  15,  1992. 


QUESTIONS  AND  ANSWERS 

J.  Barnes,  Austron,  Incorporated:  I  want  to  ask  a  quick  question  so  I  am  sure  I  understand 
what  you  are  doing  with  the  data.  You  take  a  data  link  and  take  the  first  point,  minus  the 
last  point,  plus  the  last  point,  minus  two  times  the  midpoint. 

M.  Weiss,  NIST:  Yes. 

J,  Barnes:  That  gives  you  your  estimate  of  drift. 

M.  Weiss:  Divide  by  Tau  square,  but  yes. 

J.  Barnes:  Then  you  take  that  and  calculate  a  number  to  take  from  all  of  the  data  to  get  the 
residuals. 

M.  Weiss:  Subtract  a  quadratic  based  on  that  number;  Yes. 

G.  Winkler,  USNO:  I  find  that  discussion  very  interesting.  In  fact  it  is  a  continuation  of  a 
discussion  of  drifts  which  started  about  six  or  seven  years  ago,  when  you  gave  your  paper 
about  how  not  to  measure  drift,  by  not  making  a  parabolic  fit,  for  the  phase  data;  remember 
that?  I  think  you  did  that  and  ever  since  that  time,  we  have  discussed  how  do  you  best 
measure  drift.  Before  there  was  a  question  and  it  was  questioned  whether  you  can  determine 
it  at  all,  I  believe  you  should  remember  that  whenever  we  measure  something,  we  measure  it 
against  a  hypothesis,  about  a  assumption.  You  have  in  your  various  estimates  made  various 
different  assumptions.  Each  of  these  define  — (Tape  ran  out)  — for  it’s  instability.  We  have 
to  remember  that  these  numbers  always  are  connected  with  a  assumption,  which  has  been 
made  in  the  first  place.  You  are  starting  out  with  three  point  estimates.  Why  is  it  the  best 
estimate,  or  the  optimun  estimate,  you  can  obtain,  because  it  makes  the  minimun  number  of 
assumptions. 

M.  Weiss:  I  think  you  are  right  about  the  underlying  assumptions.  In  particular,  I  think 
what  is  most  important  is  to  realize  there  is  physics  involved  and  the  reason  we  estimate  drift 
because  we  believe  that  what  is  physically  causing  the  drift  is  different  than  what  is  physically 
causing  the  random  walk.  That  they  are  two  separate  processes  and  they  should  therefore  be 
estimated  independately.  And  sometimes  random  walk  looks  an  awful  like  drift  and  there  may 
not  be  any  drift  and  it  may  be  random  walk. 

J.  Barnes:  I  see  people  like  the  idea  of  being  explicit  in  their  models.  I  think  that  is  great.  I 
think  I  like  the  comments  very  much. 
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